Differentially Expressed Genes Network Analysis

The R markdown is available from the pulldown menu for Code at the upper-right, choose “Download Rmd”, or download the Rmd from GitHub.

---

This protocol describes a network analysis workflow in Cytoscape for a set of differentially expressed genes. Points covered:

• Retrieving relevant networks from public databases
• Network functional enrichment analysis
• Integration and visualization of experimental data
• Exporting network visualizations

---

Installation

if (!requireNamespace("BiocManager", quietly = TRUE)){
  install.packages("BiocManager")
}

if(!"RCy3" %in% installed.packages()){
  BiocManager::install("RCy3")
}
library(RCy3)

Getting started

First, launch Cytoscape and keep it running whenever using RCy3. Confirm that you have everything installed and running:

cytoscapePing()
cytoscapeVersionInfo()

Prerequisites

If you haven’t already, install the STRINGapp

installApp('stringApp')

Background

Ovarian serous cystadenocarcinoma is a type of epithelial ovarian cancer which accounts for ~90% of all ovarian cancers. The data used in this protocol are from The Cancer Genome Atlas, in which multiple subtypes of serous cystadenocarcinoma were identified and characterized by mRNA expression.

We will focus on the differential gene expression between two subtypes, Mesenchymal and Immunoreactive.

For convenience, the data has already been analyzed and pre-filtered, using log fold change value and adjusted p-value.

Network Retrieval

Many public databases and multiple Cytoscape apps allow you to retrieve a network or pathway relevant to your data. For this workflow, we will use the STRING app. Some other options include:

• WikiPathways
• NDEx
• GeneMANIA

Retrieve Networks from STRING

To identify a relevant network, we will query the STRING database in two different ways:

• Query STRING protein with the list of differentially expressed genes.
• Query STRING disease for a keyword; ovarian cancer.

The two examples are split into two separate workflows below.

Example 1: STRING Protein Query Up-regulated Genes

Load the file containing the differential gene expression data and filter for up-regulated genes, TCGA-Ovarian-MesenvsImmuno_data.csv:

de.genes <- read.table("https://raw.githubusercontent.com/cytoscape/cytoscape-tutorials/gh-pages/protocols/data/TCGA-Ovarian-MesenvsImmuno_data.csv", header = TRUE, sep = ",", quote="\"", stringsAsFactors = FALSE)
de.genes.up <- de.genes[which(de.genes$logFC >= 1.8 & de.genes$FDR.adjusted.Pvalue <= 0.05),]

We will use the identifiers in the first column of this datafile to run a STRING protein query, with confidence (score) cutoff of 0.4:

string.cmd = paste('string protein query query="', paste(de.genes.up$Gene, collapse = '\n'), '" cutoff=0.4  species="Homo sapiens"', sep = "")
commandsRun(string.cmd)

The resulting network will load automatically and contains up-regulated genes recognized by STRING, and interactions between them with an evidence score of 0.4 or greater.

The networks consists of one large connected component, several smaller networks, and some unconnected nodes. We will select only the connected nodes to work with for the rest of this tutorial, by creating a subnetwork based on all edges:

createSubnetwork(edges='all', subnetwork.name='String de genes up')

Data Integration

Next we will create a visualization from the log fold changes and p-values from our TCGA dataset. Since the STRING network is a protein-protein network, it is annotated with protein identifiers (Uniprot and Ensembl protein), as well as HGNC gene symbols. Our data from TCGA has NCBI Gene identifiers (formerly Entrez), so before importing the data we are going to use the ID Mapper functionality in Cytoscape to map the network to NCBI Gene.

mapped.cols <- mapTableColumn('display name', 'Human', 'HGNC', 'Entrez Gene')

We can now integrate the same differential gene expression data with the network (node) table in Cytoscape. For importing the data we will use the following mapping:

• Key Column for Network should be Entrez Gene, which is the column we just added.
• Gene should be the key of the data(de.genes.full).

loadTableData(de.genes.full,data.key.column="Gene",table.key.column="Entrez Gene")

You will notice two new columns (logFC and FDR.adjusted.Pvalue) in the Node Table.

tail(getTableColumnNames('node'))

Visualization

Next, we will create a visualization of the imported data on the network.

setVisualStyle(style.name="default")
setNodeShapeDefault(new.shape="ELLIPSE", style.name = "default")
lockNodeDimensions(new.state="TRUE", style.name = "default")
setNodeSizeDefault(new.size="50", style.name = "default")
setNodeColorDefault(new.color="#D3D3D3", style.name = "default")
setNodeBorderWidthDefault(new.width="2", style.name = "default")
setNodeBorderColorDefault(new.color="#616060", style.name = "default")
setNodeLabelMapping(table.column="display name",style.name = "default")
setNodeFontSizeDefault(new.size="14", style.name = "default")

Before we create a mapping for node color representing the range of fold changes, we need the min and max of the logFC column:

logFC.table.up <- getTableColumns('node', 'logFC')

logFC.up.min <- min(logFC.table.up, na.rm = T)
logFC.up.max <- max(logFC.table.up, na.rm = T)
logFC.up.center <- logFC.up.min + (logFC.up.max - logFC.up.min)/2

copyVisualStyle(from.style = "default", to.style = "de genes up")
setVisualStyle(style.name="de genes up")

data.values = c(logFC.up.min, logFC.up.center, logFC.up.max)
node.colors <- c(brewer.pal(length(data.values), "YlOrRd"))
setNodeColorMapping('logFC', data.values, node.colors, style.name="de genes up")

Applying a force-directed layout, the network will now look something like this:

layoutNetwork(paste('force-directed', 
              'defaultSpringCoefficient=0.00003',
              'defaultSpringLength=50',
              'defaultNodeMass=4',
              sep=' '))

Enrichment Analysis Options

Next, we are going to perform enrichment anlaysis uing the STRING app.

STRING Enrichment

The STRING app has built-in enrichment analysis functionality, which includes enrichment for GO Process, GO Component, GO Function, InterPro, KEGG Pathways, and PFAM.

First, we will run the enrichment on the whole network, against the genome:

string.cmd = 'string retrieve enrichment allNetSpecies="Homo sapiens", background=genome  selectedNodesOnly="false"'
commandsRun(string.cmd)
string.cmd = 'string show enrichment'
commandsRun(string.cmd)

When the enrichment analysis is complete, a new tab titled STRING Enrichment will open in the Table Panel.

The STRING app includes several options for filtering and displaying the enrichment results. The features are all available at the top of the STRING Enrichment tab.

We are going to filter the table to only show GO Process:

string.cmd = 'string filter enrichment categories="GO Process", overlapCutoff = "0.5", removeOverlapping = "true"'
commandsRun(string.cmd)

Next, we will add a split donut chart to the nodes representing the top terms:

string.cmd = 'string show charts'
commandsRun(string.cmd)

STRING Protein Query: Down-regulated genes

We are going to repeat the network search, data integration, visualization and enrichment analysis for the set of down-regulated genes by using the first column of TCGA-Ovarian-MesenvsImmuno-data-down.csv:

de.genes.down <- read.table("https://cytoscape.github.io/cytoscape-tutorials/protocols/data/TCGA-Ovarian-MesenvsImmuno-data-down.csv", header = TRUE, sep = "\t", quote="\"", stringsAsFactors = FALSE)
string.cmd = paste('string protein query query="', paste(de.genes.down$Gene, collapse = '\n'), '" cutoff=0.4  species="Homo sapiens"', sep = "")
commandsRun(string.cmd)

Subnetwork

Let’s select only the connected nodes to work with for the rest of this tutorial, by creating a subnetwork based on all edges:

createSubnetwork(edges='all', subnetwork.name='String de genes down')

Data integration

Again, the identifiers in the network needs to be mapped to Entrez Gene (NCBI gene):

mapped.cols <- mapTableColumn('display name', 'Human', 'HGNC', 'Entrez Gene')

We can now import the data:

loadTableData(de.genes.full,data.key.column="Gene",table.key.column="Entrez Gene")

Visualization

Next, we can create a visualization. Note that the default style has been altered in the previous example, so we can simply switch to default to get started:

setVisualStyle(style.name="default")

The node fill color has to be redefined for down-regulated genes:

logFC.table.down <- getTableColumns('node', 'logFC')

logFC.dn.min <- min(logFC.table.down, na.rm = T)
logFC.dn.max <- max(logFC.table.down, na.rm = T)
logFC.dn.center <- logFC.dn.min + (logFC.dn.max - logFC.dn.min)/2

copyVisualStyle(from.style = "default", to.style = "de genes down")
setVisualStyle(style.name="de genes down")

data.values = c(logFC.dn.min, logFC.dn.center, logFC.dn.max)
node.colors <- c(brewer.pal(length(data.values), "Blues"))
setNodeColorMapping('logFC', data.values, node.colors, style.name="de genes down")

Apply a force-directed layout.

layoutNetwork(paste('force-directed', 
              'defaultSpringCoefficient=0.00003',
              'defaultSpringLength=50',
              'defaultNodeMass=4',
              sep=' '))

Focusing on the connected part of the network, it should look something like this:

STRING Enrichment

Now we can perform STRING Enrichment analysis on the resulting network:

string.cmd = 'string retrieve enrichment allNetSpecies="Homo sapiens", background=genome  selectedNodesOnly="false"'
commandsRun(string.cmd)
string.cmd = 'string show enrichment'
commandsRun(string.cmd)

Filter the analysis results for non-redundant GO Process terms only.

string.cmd = 'string filter enrichment categories="GO Process", overlapCutoff = "0.5", removeOverlapping = "true"'
commandsRun(string.cmd)

string.cmd = 'string show charts'
commandsRun(string.cmd)

## This part doesn't work yet
##string.cmd = "string settings chartType='Pie Chart' nTerms=3 showEnhancedLabels='FALSE' showGlassBallEffect='FALSE' defaultPalette='ColorBrewer Set1 colors' showImage='FALSE'"
##commandsRun(string.cmd)

STRING Disease Query

So far, we queried the STRING database with a set of genes we knew were differentially expressed. Next, we will query the STRING disease database to retrieve a network genes associated with ovarian cancer, which will be completely independent of our dataset.

string.cmd = 'string disease query disease="ovarian cancer" cutoff="0.95"'
commandsRun(string.cmd)

This will bring in the top 100 (default) ovarian cancer associated genes connected with a confidence score greater than 0.95. Again, lets extract out the connected nodes:

createSubnetwork(edges='all', subnetwork.name='String ovarian sub')

Data integration

Next we will import differential gene expression data from our TCGA dataset to create a visualization. Just like the previous example, we will need to do some identifier mapping to match the data to the network.

mapped.cols <- mapTableColumn("display name",'Human','HGNC','Entrez Gene')

Here we set Human as species, HGNC as Map from, and Entrez Gene as To.

We can now import the data frame with the full data (already loaded the data in Example 1 above) into the node table in Cytoscape:

loadTableData(de.genes.full, data.key.column = "Gene", table = "node", table.key.column = "Entrez Gene")

Visualization

Again, we can create a visualization:

setVisualStyle(style.name="default")

Next, we need the min and max of the logFC column:

logFC.table.ovarian <- getTableColumns('node', 'logFC')

logFC.ov.min <- min(logFC.table.ovarian, na.rm = T)
logFC.ov.max <- max(logFC.table.ovarian, na.rm = T)
logFC.ov.center <- logFC.ov.min + (logFC.ov.max - logFC.ov.min)/2

Let’s create the mapping:

copyVisualStyle(from.style = "default", to.style = "ovarian")
setVisualStyle(style.name="ovarian")

data.values = c(logFC.ov.min, logFC.ov.center, logFC.ov.max)
node.colors <- c(brewer.pal(length(data.values), "RdBu"))
setNodeColorMapping('logFC', data.values, node.colors, style.name="ovarian")

Apply a force-directed layout.

layoutNetwork(paste('force-directed', 
              'defaultSpringCoefficient=0.00003',
              'defaultSpringLength=50',
              'defaultNodeMass=4',
              sep=' '))

The TCGA found several genes that were commonly mutated in ovarian cancer, so called “cancer drivers”. We can add information about these genes to the network visualization, by changing the visual style of these nodes. Three of the most important drivers are TP53, BRCA1 and BRCA2. We will add a thicker, colored border for these genes in the network.

Select all three driver genes by:

selectNodes(c("TP53", "BRCA1", "BRCA2"), by.col = "display name")

Add a style bypass for node Border Width (5) and node Border Paint (bright pink):

setNodeBorderWidthBypass(getSelectedNodes(), 5)
setNodeBorderColorBypass(getSelectedNodes(), '#FF007F')

Exporting Networks

Cytoscape provides a number of ways to export results and visualizations:

• As an image:

exportImage('./differentially-expressed-genes', 'PDF')
exportImage('./differentially-expressed-genes', 'PNG')
exportImage('./differentially-expressed-genes', 'JPEG')
exportImage('./differentially-expressed-genes', 'SVG')
exportImage('./differentially-expressed-genes', 'PS')

• To a public repository:

exportNetworkToNDEx("user", "password", TRUE)

• As a Cytoscape JSON file:

exportNetwork('./differentially-expressed-genes', 'cyjs')